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ABSTRACT 


This study was conducted to add surface tension forces to a marker-and-cell code and 
to perform four incompressible fluid simulations in reduced-gravity. This marker- 
and-cell code has a variable grid capability with arbitrary curved boundaries and time- 
dependent acceleration fields. The surface tension logic includes a spline fit of surface 
marker particles as well as contact angle logic for straight and curved wall boundaries. 
Three types of flow motion were simulated with the improved code: impulsive settling 
in a model Centaur LH2 tank, continuous settling in a model and full-scale Centaur 
LO2 tank and mixing in a Centaur LH2 tank. The impulsive settling case confirmed 
the NASA Lewis Research Center drop tower analysis which indicated more orderly 
fluid collection flow patterns with this method providing a potential savings in settling 
propellants. In the LO2 tank, fluid collection and flow simulation into the thrust 
barrel were achieved; however, flow simulation of thin films along the bottom arbitrary 
curved boundary was inadequately modeled. The mixing simulation produced good 
results indicating both the development of the flow field and fluid interface behavior. 

The potential exists for examination of the influence of Bond number and mixer exit 
velocity on mixing time and fluid interface behavior. 
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SUMMARY 


Numerical simulations provide a valuable extension of drop tower tests to provide 
insight into fluid dynamic behavior in low-gravity environments. The purpose of this 
study was to simulate fluid behavior during settling in Centaur LH2 and LO 2 tanks 
and to examine mixing flow in the Centaur LH 2 tank using the marker-and-cell code. 

To achieve this modeling in low-g the code required modification to include surface 
tension forces. 

Initial attempts to include surface tension logic used a circular arc fit to obtain the 
radius of curvature. This approach was inadequate and necessitated going to the 
more complex spline fit for monotonic sequences. The monotonicity requirement 
of this approach was sometimes troublesome, however the method was more precise 
and the results were good. The single continuous interface has capability for 
different contact angles at each end where wall contact is made. The contact angle 
solution is usable with straight walls and concave or convex arbitrary boundary walls. 
The contact angle may be time dependent and is measured from the second or second 
last surface marker point from distances beyond the center of the adjacent boundary 
cell. The surface tension logic is basically superimposed on the surface marker 
particle logic which includes adding or deleting particles and interface smoothing. 

The surface tension force is defined as a pressure at the intersection of the interface 
and a line through cell centers but is translated to act and to be stored at the center 
of the surface cells. The method was verified with known geometric shapes where 
the radius of curvature was analytically available and was further verified with a 
test case involving fluid motion. The surface tension logic was then used on the four 
cases examined in this study. 

A second computer code improvement concerned the convergence scheme. The code 
requires extensive iterations to converge in relatively full tanks when the gradient in 
the pressure field is low in reduced-gravity fields, i.e. , coast periods rather than 
settling. At this time, the pressure field gradient in the tank is quite low. Two new 
convergence methods were considered but were found to be less efficient than the 
existing method. Therefore, the existing method was improved through an extrapolation 
technique which advances the pressure field one large step every 40 iterations. This 
method possessed stability and in all instances reduced the time for convergence. 

The code was also compiled with the FTN compiler to increase efficiency. 

Three types of flow simulation were specified under this study. Case 1 Involved 
impulsive settling with a Centaur model LH 2 tank. Cases 2 and 3 considered continuous 
settling in the Centaur LO 2 tank while Case 4 simulated mixing in a full-scale Centaur 
LH 2 tank. Case 1 was studied to provide insight into the fluid dymamics during 
impulsive settling where a fluid is set in motion by the equivalent of 12 lbs full-scale 
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thrust for a brief period after which thrust is terminated and the fluid continues to 
collect aft. A Lewis Research Center drop tower test where thrust lasted 0. 30 sec 
had indicated potential for reducing settling propellants and was simulated. The 
simulation initially progressed smoothly and flow moved down the side wall and into 
the comer. During the 0. 30 sec thrust period and until 0. 9 sec, simulation and drop 
tower results were alike. After this time in the simulation, liquid did not flow radially 
inward over the bulkhead but continued to collect in the lower comer. This was unlike 
drop tower results but not a serious deviation. Some 0. 5 sec later, flow did continue 
toward the centerline but arrived 0. 6 sec later than in the drop tower case. The flow 
velocities down the wall decreased in magnitude with time; overall behavior was as 
expected. Initial convergence difficulties in low-g made the simulation costly and it 
could not be continued until fluid motion ceased. Impulsive settling shows promise as 
a weight -saving approach to accomplish settling. Modeling can be improved if more 
rapid convergence can be achieved at low-g levels. 

In the L02 tank the behavior of settling flow in the region of the thrust barrel was of 
particular interest. Drop tower results had been only partially successful because of 
the complexity of modeling flow into the thrust barrel and also the difficulty in observing 
this flow. The variable grid capability was extensively used in simulating the thrust 
barrel. Because of cylindrical symmetry, the holes in the thrust barrel were modeled 
with slits of equal area. The model case had 30 percent liquid present to represent a 
drop tower case, while the full-scale case contained 9 percent liquid to model a full- 
scale Centaur settling. Initial flow in both simulations proceeded down the wall in a 
smooth manner. In the model Case 2, the flow proceeded up along the side of the thrust 
barrel and flow passed through the slit at a sufficient rate to indicate thrust barrel 
filling in a reasonably brief time. In both cases, however, flow along the aft bulkhead 
was difficult to model because a thin film resulted in only surface cells; these often 
lacked full or surface neighbors and could not adequately transmit the flow and pressure 
field information along this boundary. In the full-scale Case 3, an attempt to improve 
the resolution with finer cell size was inadequate as the film was significantly thinner. 
The full-scale Case 3 proceeded only to covering of the thrust barrel holes and had to 
be discontinued because of erroneous velocities along the aft bulkhead. The current 
state of the model cannot adequately simulate the flow of the thin films on the curved 
boundary. During both these simulations, surface tension calculations were satisfactory 
and convergence occurred in a reasonable number of iterations. 

Case 4 was a new area of investigation in marker-and-cell flow dynamics simulation. 

The use of a jet mixer in a low-g coast orbit was modeled for a full-scale Centaur LH 2 
tank with 65 percent liquid present. A low velocity flow was aimed toward the inter- 
face from the submerged jet. Mixing correlations predicted a time of 2000 seconds to 
achieve steady-state and total mixing while the Weber and Froude number criteria 
used predicted that interface break-up should not occur. The velocity field was 
observed to develop smoothly in the tank and progress away from the mixer and across 
the Interface. The entrained flow below the mixer was well defined. Interface 
velocities induced waves and Interface instabilities but ullage break-up did not occur. 
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Mixing times appear to be less than predicted since steady-state flow at the interface 
had occurred by 700 sec. Flow patterns indicate stratification would have been 
destroyed. Unexpected convergence problems necessitated running this case at a 
Bond number of 100 with a flat interface. Cases with lower Bond number are also 
of considerable Interest and, in view of Case 1 results at low Bond numbers, should 
be within the capability of the model. Overall results of this case were very nearly 
as expected. Since the flow field is changing only very slowly, the solution converged 
each time step in minimal iterations and the case provided a low cost simulation of 
700 seconds of mixing in reduced-gravity. 

This study indicated the potential for simulating mixing with the marker-and-cell code 
to evaluate mixing times and flow behavior at the surface. Future comparisons 
between drop tower results and model simulation would be useful in the verification 
of existing criteria for interface behavior. The surface tension capability affords 
an opportimity for analysis of behavior where this mechanism is controlling, however 
improvements in convergence methods are required for Bond numbers less than 10. 
The study also indicated a need for further code development in the treatment of thin 
films on arbitrary boundaries. 
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INTRODUCTION 


Many requirements on space missions involve the motion of propellant in a reduced 
gravity field. Drop tower tests are one method of gaining insight into this fluid behavior, 
A second method is numerical techniques such as the marker-and-cell computer code. 
These are both interim approaches until experimental space payloads are available. One 
reduced- gravity fluid behavior problem is the requirement to settle propellants to the aft 
end of a space vehicle tank prior to engine restart. During space coast in earth orbit, 
propellants collect forward due to atmospheric drag or retromaneuvers and present a 
settling requirement. Propellants aboard Centaur are settled with the application of 
a settling thrust from hydrogen peroxide motors. The pattern of fluid motion during 
settling is of interest to define the required settling times. The operation of a jet 
pump in a cryogenic propellant tank in reduced gravity is a second fluid problem. A 
fluid jet has been recommended for mixing the propellant to reduce pressure rise due 
to stratification at the vapor interface. The extension of one-g mixer data to low 
gravity conditions has raised some doubts however, and until very recently no low'-g 
drop tower data was available in this area of reduced-gravity research. 

The marker-and-cell technique provides a convenient analytical tool to simulate low-g 
fluid behavior and to define both flow and pressure fields as well as the time require- 
ments for the motion. Visual presentations of the fluid location and the velocity field 
are provided as well as detailed printouts on the velocity and pressure field. In very 
low-g fields, the surface tension forces dominate the acceleration forces. Therefore, 
it was important to add logic to the program which solves for these forces on the inter- 
face, particularly for simulation of impulsive settling or of mixing where acceleration 
forces are much lower than ejq>erienced in normal continuous settling. 

In each of the cases considered in this study, drop tower data were available for model 
tank configurations. This provides for verification of a numerical method which can 
then be used for analysis of full-scale simulations and for extension of simulations 
beyond the time-limits of drop tower experiments. The logic and capabilities of the 
code in its original form and the logic for the surface tension modifications are both 
outlined in this report, as well as simulation results. 


PRECEDING PA.GE BLANK NOT 




COMPUTER CODE DEVELOPMENT 


The computer code development task in this study had two objectives; the modification 
of a marker-and-cell code to include surface tension and the improvement of convergence 
techniques in the code when simulating very low gravity fields. An existing marker-and- 
cell code, ERIE (Ref. 1-4) was used as a starting point for this modification. The initial 
code contained a variable-grid capability and was in overlay structure to reduce core 
storage. 

The code Is used to solve problems in incompressible fluid d3mamics using an approach 
known as the Simplified Marker-and-Cell (SMAC) technique which was originated at Los 
Alamos Laboratory (Ref. 5). The cells provide a grid mesh for problem solution and 
liquid accountability with a specification of velocity and pressure fields; the particles 
indicate the fluid location. The newly improved code has the capability for modeling 
surface tension pressure, arbitrary boundaries, and variable gravity fields. 

In order to acquaint the reader with the basic approach to the computer logic and to 
provide stand alone capability for this document, a few pages from Reference 1 are 
summarized to present further insight into the marker-and-cell approach prior to 
describing modifications current to this contract. For the reader interested in greater 
detail on the program logic, a current program documentation is provided in Reference 6. 

The SMAC method is a numerical finite difference technique for solution of the Navier- 
Stokes equations of motion for incompressible viscous fluid flow. The method solves 
the complete Navier- Stokes equations givir^ a time-dependent solution. The method is 
applicable to confined flow or free surface flow. The SMAC method overlays the fluid 
with a two-dimensional rectangular grid mesh in either rectangular or cylindrical 
coordinates. Horizontal velocities are defined at the midpoints of each side of each 
rectangular element while vertical velocities are defined at the midpoints of the top 
and bottom sides. Finite velocities are defined each time step only for cells which 
are full or surface; the latter being those next to an empty cell. Various designations 
are given to other cells in the grid depending on their location and function as shown 
in Figure 1. Velocities satisfy tangential and normal stresses on surface cells. 

A pressure field is defined throughout the fluid with the pressure defined at the center 
of each full cell. The pressures on surface cells are interpolated to satisfy the normal 
stress condition plus a surface tension pressure force. The calculation proceeds with 
a calculation over all full and surface cells during which velocities are set independent 
of pressure and the divergence requirement. These velocities assure that vorticity is 
satisfied at each cell comer. The iteration procedure in the SMAC code then modifies 
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ARB 

BND 

BOR 

BOX 

COR 

BMP 

EOC 

EXT 

PUL 

LEP 


OB 


OK 


RIG 


SNC 

SUR 

TOP 


An OB cell which has a fluid particle within t - DR of the arbitrary boundary, where c is the boundary 
sensing parameter and is usually set equal to 0.25f 


A border cell which is neither COR or EXT , 

A cell bordering the interior cells- A border cell is not an interior cell and cannot contain fluid . 

A COR cell containing a segment of an arbitrary boundary which has its midpoint and angle being 
stored in the OB cell Just below it. 

A cell which has a line segment of the arbitrary boundary passing through it. however, fluid area to 
total cell area fraction is less than 0-25. The fluid area is to the left of the line segment. Also, any 
cell just outside an OB cell is a COR rell- 


The cell is empty (contains no fluid 
particles). 

A coll which is either EMP or COR. 

Any cell outside a COR cell. 

A cell which contains fluid and has 
no empty neighbor. 

A COR cell containing a segment of 
an arbitrary boundary which has its 
midpoint and angle being stored in 
the OB cell Just to the left of it. 

A cell which has a line segment of 
the arbitrary boundary associated 
with it (see COR). 

An ARB cell that does not contain 
the Intersection of a free surface 
and the boundary. 

A COR cell containing a segment of 
an arbitrary boundary which has its 
midpoint and angle being stored in 
the OB cell Just to the right of it. 

A cell which is also flagged as FUL 
or SUR and not COR. 



A cell which contains fluid and has at 
least one neighbor with an empty 
neighbor not COR or EXT . 


1. All ARB cells are also flagged as 0B. 


A COR cell containing a segment of 
an arbitrary boundary which has its 
midpoint and angle being stored in 
the OB cell just above it. 


2, All T0P, B0T, RIG and LEF cells are also flagged 
as C0R cells. 

3, All the above SUR and FUL cells are also flagged SNC, 


Figure i. Variable Grid Network Illustrating Typical Cell Flags 
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these velocities to minimize divergence over all cells by adjusting the pressure in each 
full cell by an incremental amount. The iteration is complete when a specified 
convergence criteria is satisfied based on a measure of the pressure change which 
is proportional to divergence. 

When convergence is satisfied, the velocities calculated above are used to move an 
array of massless particles defining both the surface and the fluid location. It is 
these marker particles which give the method its name, marker and cell. The surface 
is defined not only by the above criteria for surface cells but also with an array of 
surface particles which are moved in a similar manner to marker particles. These 
surface particles are also used to determine the magnitude of the surface pressure. 

The graphic output of these particles each time-step and a similar printout of cell- 
centered velocity vectors are the descriptive output of this method. This is in addition 
to velocity, pressure, divergence data, and cell flag data which are printed out for 
each cell. 

2 . 1 EQUATIONS OF MOTION 


The computer code ERIE solves the Navier-Stokes equations and satisfies the equation 
of continuity, i.e. , divergence for full fluid cells is zero. The basic differential 
equations are presented below 


^ J_ 5r u 5uv 

Bt a 5r 5z 

r 


90 

9r 


+ Sr + u 


9z 




-s i -s O' >2 

3v 1 Sr nv 3v 

^ 

9t a 9r 9z 

r 


90 . u 9 

— +g (t) 

9z Of 9r 

r 




Of 9r 9z 


( 1 ) 

( 2 ) 

(3) 


The pressure 0 is normalized with fluid density and has xmits L^T"'^. The above 
equations are applicable to rectangular coordinates when o; = 0 and apply to cylindrical 
coordinates when ck- = 1. The SMAC code utilized ZIP-differencing techniques to solve 
the above equations using variables identified at specific cell locations. The modification 
of that code to variable grid, VGSMAC, introduced pecularities to the differencing 
techniques. 


A cell is flagged as a surface (SUR) cell when it contains fluid marker particles and 
it has at least one adjacent neighboring cell which is flagged empty. On free surfaces 
the tangential stress condition is 


9 u ^ 9 V 
■9 z 9 r 


(4) 
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In addition, the normal stress condition is 

0N = % (applied) - vg) (5) 

The applied part of the pressure is specified according to the requirement of the 
problem while the viscous part assures that there is otherwise no net flux of normal 
momentum through the surface. It is important that the normal stress condition be 
placed on the free surface rather than at the center of the surface cell. 

Tentative velocities are calculated for all Ml and surface cells without the pressure 
contribution. These tilde velocities are determined before the final iteration scheme. 
Then, the pressure in each full cell is modified to minimize the summation of a 
normalized divergence for all full cells. Thus, tentative velocities are modified by 
Equations 6 and 7 in the iteration scheme 



The viscosity coefficient u used in Equations 1 and 2 is the sum of a kinematic molecular 
viscosity and a turbulent viscosity 


V 


u 


molecular 


+ 


u 

turb 


( 8 ) 


The molecular viscosity is an input quantity and is a fluid property. The turbulent 
viscosity coefficient is calculated internally in the program as indicated below 



This expression for turbulent viscosity is of the form predicted by both Prandtl’s mixing- 
length theory and Taylor's vorticity transport theory. A turbulent viscosity is calculated 
in a cell containing fluid when at least two of its adjacent neighboring cells also contain 
fluid. The criteria of requiring fluid in adjacent fluid cells is needed so that Sv/Sr and 
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ou/Bt; can be calculated. For the reorientation flow cases examined in this study 
(Section 3) the magnitude of the local turbulent viscosity coefficient was at least an 
order of magnitude greater than the molecular viscosity coefficient for most of the 
problem duration. This indicates that the viscosity coefficient used in Equations 1 and 
2 is mainly a result of the turbulent viscosity coefficient. 

2. 2 ERIE COMPUTER CODE 


This section describes in detail the calculational cycle in the framework of ERIE. 
The code ERIE contains a number of features that make it a useful analytical tool for 
fluid dynamic studies. In particular, 

a. The code is written in FORTRAN IV for the CYBER 172 system in overlay 

structure (Figure 2) and interfaces with the SC-4020 plotter to provide particle 
and velocity vector plots. 


b. The code can be used with either cylindrical or rectangular geometry. 
Axisymmetric flow can be computed in a cylindrical geometry mode. 


c. Grid mesh si7.e may be varied for different areas of the model with a variable 
grid. The size of the computing mesh can be changed from problem to problem, 

including both the number and 
size of cells and the number of 
particles. The setup allows 
for different initial conditions 
and particle resolution in 
specified regions of the mesh. 

A scaling feature provides for 
model and full-scale simulation 
by changing only one variable. 



Figure 2. Basic Logic Flow of Overlay 
for Program ERIE 


d. Various boundary conditions 
are available. Both curved 
(arbitrary) and straight wall 
boundaries can be used to 
represent various shapes. 

e. This code can calculate both 
free-surface and confined flows. 
Surface conditions may include 
surface tension forces. 

f. Gravitational effects in 
rectangular coordinates may 
be included in any orientation 
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and may be time-dependent. In cylindrical coordinates, only the axial gravity 
force is appropriate due to symmetry; it may be also time-dependent. 

2. 2. 1 OVERLAY STEUCTURE - Program ERIE makes maximum usage of overlay 
structure in that the time-consuming pressure iterations are done in a link (3,0) 
requiring low core storage (Figure 2). The surface tension logic added during this 
study significantly expanded link (2,2); however that link is not used but a small 
fraction of the time. The velocity and marker particle plots make this a most 
descriptive computational tool. The amount of data generated cannot be adequately 
assessed without reviewing the plots which are generated. 

2.2.2 VARIABLE GRID MESH AND CELL FLAGS - The variable grid concept 
affords a method of analyzing particular sections of the problem in finer detail. 

Each cell is rectangular as it is in the fixed grid model from which the variable 
grid program was derived. Each interior cell may have either one or two neighbors 
in each of the four directions, left, top, right or bottom, i.e. , a maximum total of 
eight neighbors. The primary cell variables in the ERIE code are velocities u and v 
on right and top sides, respectively; the pseudo -pres sure '[f which is cell-centered 
except for boundary cells where it is located at the segment midpoint; and the velocity 
divergence, D, which is satisfied for full and surface cells and is the convergence 
criteria on full cells. The calculation of all these variables is made more complex 
by the unequal lengths of common cell sides. Special provisions were taken to 
assure cell variables were interpolated at the correct locations outside the cell 
being determined. This program computes the above variables throughout a grid 
mesh in which cells are flagged with one or more flags indicating their fluid state 

or relation to the boundaries defining the problem. These cell flags were defined in 
Figure 1 where an example is presented of a variable grid configuration with 
appropriate cell flags indicated. Only the most significant flag is shown on a cell; 
a cell may be flagged 0B, ARB, FUL, 0K, and SNC simultaneously. 

In addition to the mesh of Eulerian cells, the ERIE code employs a set of massless 
marker particles which provide a visual representation of the fluid. An essential 
purpose of these marker particles is to define the position of full and surface cells so 
that the configuration of the surface can be sensed. Beyond this, the marker particles 
do not enter into the calculation, but are merely embedded in the fluid and are carried 
along by it. The particles are inserted at a density per given problem dimension. 

Since particles primarily impact the flagging procedure, it is often desirable to input 
denser particle arrays where fluid motion is most rapid or is most interesting. Each 
cycle the marker particles are moved with an area weighted average of the four nearest 
u's and of the four nearest v’s. 

2.2.3 BOUNDARY CONDITIONS - In ERIE both straight line and curved wall boundaries 
can be used. The ERIE version of the program uses techniques developed by Vieeelli 
for curved-wall arbitrary boundaries (Ref. 7). A series of points are input to specify 
the curved boundary. These points are a series of coordinates which do not necessarily 
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lie on the grid-mesh lines. To define the arbitrary boundary segments on the grid 
mesh, successive points are connected by straight lines and their intersections with 
the grid mesh are determined. Only one straight line per cell is permitted with 
intermediate points being dropped. 

The following relaxation equation is used to compute the pressure in full arbitrary 
boimdary cells 




k+1 

N 



RELAX 

X 



n+l\k 
M J 



{ 10 ) 


In this equation n is the normal at the midpoint of the boundary segment associated with 
cell N and (Vm’’ ' is the liquid velocity at the midpoint of the segment computed with 
an area-weighted interpolation formula. Clearly iterates and 

must be recomputed each time the pressures and velocities are adjusted. The relaxa- 
tion parameter and minimum mesh dimension are RELAX and X, respectively. Note 
that pressure is not adjusted to satisfy divergence or net flux out of the cell but is 
adjusted proportional to the velocity directed perpendicular to the normal vector at 
the segment midpoint. If liquid is flowing across the boundary the pressure is increased 
until the outflow stops. If liquid is tending to separate from the boundary the pressure 
is decreased until the liquid flows tangent to the boundary. 


Where free surfaces intersect curved wall boundaries, special procedures are required. 
The pressure at the intersection point should equal the ambient pressure. However the 
pressure is defined only on the Eulerian net and it is sometimes not possible to zero the 
flux at the boundary consistent with vanishing divergence without introducing a pressure. 
In most circumstances the liquid in the cell will be part of a much larger mass. When 
there are one or two liquid neighbors, the velocity components at the sides in contact 
with the liquid are preserved, and those at the open and boundary sides adjusted to 
make the velocity tangent at the boundary consistent with vanishing divergence. In the 
case of one liquid neighbor, the velocities at the opposite cell sides are assumed equal, 
and the component with both sides open or boundary is adjusted. If the velocity at the 
boundary is initially directed away from the boundary, nothing need be done. The 
remaining possibility is that there are no liquid neighbors, as happens when a small 
isolated element strikes the boundary. In this case we set the component of the particle 
velocity noimal to the boundary equal to zero and preserve the tangential component. 

If a gravitational force is present we accelerate the particle velocities by the component 
of the gravitational vector tangent to the boundary. This is a free slip condition. 


2, 2. 4 TIME-STEP - The procedure for calculation of maximum time-step uses two stabil- 
ity and two accuracy criteria. The length dimensions (DRSTEP, DZSTEP) are selected 
equal to the average of the cell mesh sizes, which may exceed the accuracy criteria of 
the smallest cells. Therefore, care must be taken in selecting these variables which 
determine the maximum allowable time-step. The calculated time-step is used except 
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that when a print interval time is exceeded, the time- step is dropped to the print time. 
To avoid an extremely small step in this latter case, the criteria is relaxed so that a 
step to a print cycle may exceed the calculated time-step by forty percent. 

The first stability criteria exists to satisfy conditions imposed by the differencing 
technique 


6ts 


(DRSTEP)^ (DZSTEP)^ 

4 p [(DRSTEP)^ + (DZSTEP )^] 


( 11 ) 


The second stability criteria exists to satisfy *'Courant" condition and uses the wave 
speed and fluid depth 


2 * DRSTEP *DZSTEP 

^ — — ( 12 ) 

(DRSTEP +DZSTEP) c 

where c is the wave speed. This stability criteria was not a limiting value during this 
study. With active surface tension, c = ^gh + 27r(^3 • DRSTEP.’ 


The third and fourth restrictions are related to accuracy so that a particle does not 
move across an entire cell in a time-step. These accuracy criteria were the limiting 
criteria during this study 


DRSTEP 

6t 

2u 

max 

DZSTEP 

— 

2v 

max 


(13) 

(H) 


To insure that particles do not move across an arbitrary boundary during a time step, 
the particle sensing parameters are included 

8t ^ DRSTEP X 2 X €/u (15) 

max 

6t ^ DZSTEP X 2 X e/v (16) 

max 

where e is the boundary sensing parameter. This latter criteria is only used when the 
fluid is going toward the boundary (not parallel to or away from the boundary). 

2. 2. 5 CONVERGENCE CRITERIA - To reach an acceptable solution, the iteration 
scheme continues to adjust the pressure throughout the fluid, i. e. , the FUL and OK 
cells, until the convergence criteria have been satisfied. Recall the tilde (temporary) 
velocities are adjusted by the pressure gradient to obtain the final velocity distribution 
which is tested for divergence. Convergence criteria are tested against the sum of the 
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normalized divergence over all cells rather than the divergence of a single cell which 
would be a more stringent requirement. The tilde velocities are modified into final 
velocities for the cycle in a way that will preserve the vorticities that have been 
correctly implanted into the fluid, hi the case of non- OK cells this will now bring 
the divergence to zero, while in the case of OK cells it will cause the fluid to flow 
parallel to the arbitrary boundary. A delta pressure is determined for each FUL 
and OK cell for each iteration. New velocities are then computed for that cell as a 
function of the newly defined pressure before preceding to the next cell. When all 
cells have been treated in this fashion, a test is made to see if this pressure field 
resulted in convergence. If not, the procedure is repeated for another iteration. 

At each iteration, k, the change in pressure for FUL and OK cells is a function of the 
quantity D (D = divergence for FUL cells) which is computed for all FUL and OK cells 
according to the following expressions, 

D Calculation for FUL Cells; 


D 


.k+1 

N 


“■n l 


k k+1 

flip u - RIP u^^/ 
N N NL NL 


&z 


N 


k k+i 
V - v 
N NB 


(17) 


where RIPn and RI^NL radial distances to cell edges. 
D Calculation for 0K Cells: 

DN0RX^ * UK^ H- DN0RY^ * VK^ 


D 


N 


MAX (6r^, 6z^) 


(18) 


where UKj^ and VK are the u and v components of the velocity at the midpoint of the 
boundary segment associated with cell N. 

The pressure iteration proceeds through the cells in the order in which they were input. 
Only one matrix is used for i[t, so that the program uses the latest values (iteration 
number k+1) for velocities whose index is smaller than that of cell N and old values 
(iteration number k) for those velocities whose index is greater than N. 

The pressure of cell N at the end of k+1 iterations is then defined by 


k+1 k k+1 

% ’^N " % 


(19) 


where 5 N ® relaxation parameter defined by 


ALP {[MIN (6r^, 6z^)]^} 


^N 


N' 


6t 


( 20 ) 
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where ALP is an input variable. The convergence test is 


E 

N=1,NCELL 

N=1,NCELL 


TEMP2 


ll/2 


(^^-TEMP)^ 


< EPS 


( 21 ) 


where TEMP = §x for FUL cells has been satisfied for all FUL cells. Here EPS 
is usually on the order of 8 x lO”"^ x MIN (6rNl/2, 8 zn 1/2). (Ref. 7). 


On the relaxation of the gravitational field in impulsive-g settling, convergence problems 
were encountered. The pressure field is the key variable in convergence; a field which 
ranged from 0 to 1000 was decreased to the order of 0 to 10 when the liquid head was 
renaoved. Several convergence techniques were attempted to reduce the number of 
iterations. One scheme involved convergence of single cells rather than overall 
convergence. Another scheme involved a sophisticated technique encompassing a 
steepest ascent gradient method. The latter proved to be less efficient in this situation 
than the current approach. 


To reduce the iterations, an approach valid in the solution of certain classes of simultan- 
eous equations was employed. Three consecutive iterations are used to establish a 
trend for extrapolation of each pressure value to a new value. The equations used are 


(k) _ (k-1) 

^new_^(k-l) ^ J i 
i i 1-// 

in which fi is defined by 



(22) 


(23) 


In the preceding equations, the superscripts k-1, k, and k+1 identify the three previous 
Iterations, k+1 being the most recent. The variable p is constrained to be = 0. 99 to 
prevent overextrapolation. 

Currently this procedure is not activated until 60 iterations have passed. It is neces- 
sary to wait a sizeable number of iterations to assure that transients have disappeared 
and convergence is underway. Additionally, extrapolation is repeated every 40 
iterations after the initial 60 if required. The above technique has been checked out on 
the impulsive-g model case and in all instances convergence was attained. This 
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procedure resulted in improved convergence in comparisons made with and without it. 

In no case were the results slower. It was made a permanent part of the program and 
used in all four cases of this study. Frequently, convergence was reached in less than 
sixty iterations on Cases 2 and 3. For Case 4, a typical number of iterations was 12, 
the minimum permitted for this case. It has been recognized that less than 12 iterations 
may lead to false solutions. 

2.2.6 SURFACE TENSION FORCE - A major task in this study was the improvement 
in the surface tension logic used in computer code ERIE. In investigating fluid dynamic 
behavior at reduced gravity, the surface tension forces become important at Bond 
numbers below 10. Below Bond numbers of 1, surface tension forces dominate accelera- 
tion forces. Since fluid behavior during the coast phase of the Centaur vehicle is of 
interest when Bond numbers are 10 and below, it is necessary to consider surface 
tension forces. 

In this study a method suggested by Daly (Reference 8, 9) was used for adding 
surface tension forces to the marker-and-cell code. The surface tension force is 
applied as a pressure across the interface acting at the intersection of the interface 
with lines through surface cell centers. In the following paragraphs, various details 
of the logic will be discussed. 

The surface tension pressure is defined by the equation 


'^ST ~ ^ (l/^l 1/^2^ 


(24) 


where jSis a kinematic surface tension coefficient (L^/T^) and Ry is a radius of curvature 
in the two dimensional plane of the paper (a transverse force) and R 2 is a radius of 
curvature (a circumferential force) in a plane peipendicular to the former plane and 
passing through the interface point perpendicular to the tangent to the curve in plane 
one. In plane coordinate systems, R 2 is not considered. 

This method makes use of the surface marker particles which define the surface 
location. These are moved each time step with the velocities which have been determ- 
ined for each full or surface cell. After movement of the particles, the pressure field 
along the interface is determined from the curvature of the interface defined by these 
particles. Particles are checked for overall spacing to assure that movement has not 
resulted in large gaps or particles moving too close to one another. If this occurs, 
particles are added or deleted as the situation requires. 

The next step in the process is the relocation of the end points of the surface line. 

These points will have been moved earlier with the local velocities. Since no-slip 
boundaries are used, the particles will move parallel with the straightline boundary 



with the velocity of the adjacent cell. For arbitrary boundaries, a similar parallel 
flow is assured. However, the location of each end boundary point is still relocated 
each step to satisfy contact angles which are input as a specification j contact angle == 
consti + const 2 (time)®®^®*'3 for each end. Boundaries are input as a series of points 
defining straight lines; as few as two may be input to define a line for intersection of 
each end of the interface. Although zero contact angle fluids may be used, the 
macroscopic angle from the center of the adjacent cell is larger and this is the angle 
which is required. The proximity to the boundary of the second and second to last 
surface marker particle is limited by an input specification of a peipendicular distance 
to the wall, (Figure 3). If the particle is closer than this distance, the particle is 

deleted and the next particle is tested imtil one is 
found which satisfies this minimum distance 
requirement. The contact angle is computed 
between this point and the appropriate boundary 
segment and a point defined on this boundary 
segment which satisfies the angle requirement. 
Special logic was required for convex boundaries 
since perpendiculars and some specified contact 
angles are not always defined. In these cases, 
the nearest point on the boundary is located and 
used to satisfy the boundary point for that time 
step. 

After the points have been defined, they are divided 
into X and y sequences, since spline fits require 
monotonicity. Spline fits also require five points in 
the sequence. Therefore, if five sequential points fail to satisfy the monotonicity test, 
smoothing is applied. If no improvement occurs, a circular arc fit is used. The sequence 
switches from r to z when the slope of five consecutive segments exceeds 1. 5 or from z 
to r when the slope is less than 0.666 for five consecutive particles 


n-4 



Figure 3. Contact Angle Definition 
From Interface Particle (n-1) to 
Particle (n) 


and 



> 1.5, 


i — k, k 1, k + 2, k + 3, k +4, 


\ \-l 

\ ' Vi 


<1.5 


(25) 


(26) 


A sequence of points with a tight radius of curvature are also smoothed to avoid 
unrealistic pressures due to local deviations. The spline fit, which uses five points, 
also requires monotonicity of two points into the next sequence, with program stops 
resulting if this is not satisfied. An exchange of two surface particles sometimes 
satisfied this requirement. 
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After the sequences are defined, the sequences are solved as a series of cubics with 
a spline fit interpolation scheme. It is the coefficients of the matrices which result 
that are used to define first and second derivatives at the r, z locations of the inter- 
section of the interface segments and lines through surface cell centers. The radius 
of curvature are defined by the following equations for the radial sequence 


Rl - [{z')2 + z” (27) 

Rg -r[(z')2 + l]^'^y jz' 1 (28) 

and for the axial sequence 

Rl = [(r')^ + ^ r" (29) 

R 2 = r [ (z’ )2 + Ir' 1 (30) 


With the radius of curvature defined, the equation for surface tension pressure is 
solved. In the Marker-and-Cell technique, pressures are always defined at cell 
centers with the exception of arbitrary boundary cells. Thus the pressure must be 
translated to the center of the surface cell based on the following equation, as was 
illustrated by Figure 4 taken from Reference 10 

%={l-n)^F + ('^ST + 't'NS)T) (31) 

in which 

tfg = the pressure at the center of the SUR cell 

i|tp - the adjacent FUL cell pressure (e. g. point C in Figure 4) 

it'gT = surface tension pressure at the closest 
intersection point as described above 

'I'NS ~ normal stress pressure 


The value of T)is confined to 0.667 ^ T) = 2.0 to insure 
that only pressures from nearly full cells are used. 
Equation 31 is used so that when the regular formulas 
are used for computing velocities in FUL cells adjacent 
to SUR cells, the desired pressure is simulated at the 
surface rather than at the center of the SUR cell. 
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Figure 4. Surface Tension 
Pressure Solution 


19 




2. 2.7 SURFACE TENSION VERIFICATION - A series of runs were made to verify 
the surface tension logic. A circular interface was input and the radius of curvature 
was evaluated at the intersections of cell centers to verify the calculation results, A 
similar input was made for a sine curve so that the known values could be compared 
with the program calculated values. The two surface shapes indicated appropriate 
values were being calculated for the surface tension pressure. 

A typical settling fluid configuration was input to evaluate the method of wall contact 
angle simulation for both straight line and curved boundaries. The fluid was placed 
forward in the tank to simulate conditions prior to settling for engine restart. The 
fluid is Initially in contact with the cylindrical wall and later with the tank aft bulkhead. 
Results of this simulation are shown in Figure 5. This particular sequence is shown 
to indicate the transfer of the end point data from one boundary to the other. Until 
0. 2435 seconds, the interface remained attached to the vertical straight line boundary 
(cylinder wall) because the end point data only included two points on that line. In a 
problem restart, interface end point boundary information for the arbitrary boundary 
(bottom bulkhead) was input. The interface then transferred automatically from the 
side to the bottom boundary and the problem continued smoothly. 


20 










3 

FLUID DYNAMIC SIMULATIONS 


The marker-and-cell program ERIE with modifications for surface tension (Section 
2 . 2. 6) was used for the simulation of both model and full-scale fluid motion, A total 
of four cases were considered, two model settling cases modeling the Centaur fuel and 
oxidizer tanks, a full-scale settling case for the Centaur oxidizer tank, and a mixing 
simulation in a full-scale Centaur fuel tank. The former two cases duplicated similar 
investigations in drop tower tests at the NASA Lewis Research Center. The fuel tank 
model case was unique in that the potential for impulsive settling was investigated. 
Details of these cases are summarized in Table I. 


TABLE I. FLUID AND PROPERTY DATA FOR FOUR MODEL CASES 



— 

Radius 

cm 

— 

Fluid 

% Liquid 

Bond 

Number 

Initial 

Accelera- 

tion 

cm/sec^ 

V 

cmVsec 

(8 

cm^/sec^ 

Case 1 

5.5 

Ethanol 

65 

15 

-73.5^ 

0.01520 

28.3 

Case 2 

7.5 

FC-782 

25 

10 

-60.8 

0, 00477 

7.67 

Case 3 

152.4 

LO 2 

9 

15 

-1.10 

0. 00167 

11.7 

Case 4 

152.4 

LH2 

65 

100 

+0. 116 

0. 00192 

27.0 


Note 1. Acceleration set to +0.001 cm/sec^ at 0.30 sec after impulsive settling. 
Note 2. A fluorocarbon solvent registered by Minnesota Mining Mfg. Co. 


The initial fluid velocities for all cases were zero. The interface configuration was 
defined by the initial Bond number stated in Table I with the exception of Case 4 where 
a flat interface was selected. Interface shapes for cylinders were determined from the 
data by Hastings (Ref. 11) and the interface shapes in the oxidizer tank were from 
Concus (Ref. 12) who presents data for ellipsoids. The settling accelerations resulted 
in Bond numbers of 80, 450, and 1800 for Cases 1 through 3, respectively. For Case 4, 
the initial Bond number of 100 was held constant for the duration of the run. 

Other problem specifications common to all cases included free-slip boundaries and a 
turbulent viscosity coefficient of 0, 05 (Ref. 3). All cases were modeled with only the 
right half of the tank represented, assuming S5rmmetry about the centerline; the problems 
were run with the cylindrical coordinate option. Although the Centaur fuel tank has 
baffles, neither fuel tank simulation here considered baffles. The drop tower model 
did not use baffles. Baffles are within the capability of marker-and-cell techniques 
and were considered in Reference 1. In going from a model to a full-scale tank, the 
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same cell configuration can be used and the change in dimensions covered by a scaling 
factor. This was done for both full-scale cases where tanks of 7. 0 and 6. 696 cm were 
scaled by factors of 21. 7714 and 22. 7584, respectively. A scaling factor was also used 
to reduce Case 1 from 7. 0 to 5. 5 cm. The number of cells and particles used will be 
discussed with the individual cases, 

3. 1 IMPULSIVE-G SETTLING - MODEL CASE 1 

Impulsive settling consists of applying thrust to the vehicle for a brief period of time to 
set the fluid in motion, then stepping down or ceasing thrust until such time as the fluid 
is nearly collected. When that occurs, settling thrust is again applied prior to engine 
start or engine start may occur with no further settling thrust. The NASA Lewis 
Research Center performed a drop tower test to demonstrate the potential for this 
approach to propellant settling. The results of their test are illustrated schematically 
in Figure 6 where the interface shape is shown at various times during the drop. Note 
that the liquid has not reached the bottom bulkhead when thrust ceases at 0. 30 sec, yet 
the fluid continues to collect until approximately 3 seconds. Beyond this period the fluid 
returns to a static interface configuration. 

It was the object of this study to model the above described behavior with the ERIE 
code. Surface tension effects dominate in this problem after the thrust is terminated, 
presenting a valid test of surface tension forces, A Centaur model fuel tank of the 
dimensions used in the drop tower test, similarly without baffles, was selected 
(Figure 7). The configuration utilized 400 cells with an initial array of 2272 particles. 
Final particle array was 2728 particles. The flow initially moved down the side wall 
at a rate similar to the drop tower test. Marker particle and velocity vector plots 
are shown in Figures 8 and 9. The grid configuration is shown in Figure 9A. This 
was the first simulation with surface tension active. A contact angle of 20° from the 
center of the cell nearest to the wall was selected to represent this wall phenomena. 

This was derived from a consideration of the macroscopic angle which results for 
settling Bond numbers with zero contact angle fluids. 

Flow simulation went smoothly down the sidewall, the flow reaching the comer at 0. 5 
seconds, not appreciably different from model results. A filling-in occurred in the comer 
area and flow simulation continued to follow the model at 0. 86 seconds elapsed. Veloci- 
ties in a radial direction along the bulkhead were weaker than anticipated from this time 
onward. After 0.9 seconds, velocities in the lower comer region took on an eddy motion 
and curled back over the top to fill in fluid in this region rather than proceeding across 
the bulkhead. This phenomenon continued for nearly one-half second and is not explain- 
able. During this period the simulation falls behind the drop tower model. At approxi- 
mately 1. 5 seconds the flow again continued toward the centerline. The simulation was 
terminated at 1.91 seconds when flow had nearly reached the centerline. 

One measure of the progress of this simulation was a slowing down of velocities along 
the wall due to turbulent dissipation. This phenomena started at 0. 3 sec when thrust 
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Figure 6. Drop Tower Results for Impulsive-Settling With Model Centaur LH 2 Tank 
Case 1 - 0, 30 Seconds Thrusters On 
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Figure 7. D-IT Tank Configurations With Full-Scale and Model 
Dimensions 





A. T=O.OSEC B. T=0.15SEC C. T= 0.30 SEC D. T=0.62SEC E. T = 0.625SEC F. T = 0,80SEC 


G. T = 0.95 SEC H. T= 1. 15 SEC I. T= 1.34 SEC J. T= 1. 55 SEC K. T= 1.76 SEC L. T= 1.88 SEC 



Figure 8. Marker Particle Plots for Model Simulation of Impulsive-G 
Settling, Case 1 
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A. T= 0.0 SEC B. T= 0.15 SEC C. T = 0, 30 SEC D. T= 0.52 SEC E. T= 0.625SEC F. T= 0.80 SEC 



G. T = 0.95SEC H. T=1.15SEC I. T=1.34SEC J. T=1.55SEC K. T=1.76SEC L.T=1,81SEC 

Figure 9. Velocity Vector Plots for Model Simulation of Impulsive-g 
Settling, Case 1 
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terminated and was tracked tor several sidewall cells where the thin fi lm existed for 
tlie duration of the run. The results are shown in Figure 10 where velocities down 
the wall are plotted versus simulation time; cell distances above the bottom comer 
are identified. It was anticipated the problem could be continued until such time as 
surface tension effects caused the flow pattern to stagnate and reverse to an upward 
direction as occurred in the drop tower model (Figure 6), The flow stagnation 
described above and convergence problems prevented the evaluation of this behavior. 

Difficulties with convergence were encountered during much of this run beyond 0. 30 
seconds. This problem was discussed in Section 2 where a temporary solution to this 
difficulty was presented. 

3. 2 CENTAUR LO2 TANK SETTLING - MODEL CASE 2 

This simulation was a settling study of the Centaur oxidizer tank. It was a case pattern- 
ed after drop tower work by the NASA Lewis Research Center. The Centaur full-scale 
oxidizer tank is shown in Figure 11. The thrust barrel was a major item in establishing 
this model. The basic structure and openings in the full-scale thrust barrel are 
shown in Figure 12C. In this 1/20 scale drop tower modeling, the holes were simulated 



Figure 10. Velocity Decay Down Sidewall in Impulsive Settling, Case 1 
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A. Model Simulation 
Dimensions - Case 2 


B. Full“Scale Simulation 
Dimensions - Case 3 


C. Full-Scale 
Dimensions 


Figure 12. Centaur D-IT LO 2 Tank Thrust Barrel and Cover Details 
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with only two slits which, in reality, are annuli in a three-dimensional model due to 
axial symmetry. A top slit represents 8 of the 0. 5 cm dia holes while a side slip 
represents 20 of the 0. 3 cm dia holes. In the model, the tank radius is 7. 5 cm and 
the height 11. 0 cm; the eccentricity of the tank is 0. 68. The radius of the thrust 
barrel. Figure 12A, is 3.25 cm, the height is 2.0 cm. The side holes are represented 
with a 0. 25 cm slit and the upper holes with a 0. 25 cm slit. 

The tank was modeled with 239 cells with 812 particles indicating the 25 percent liquid 
residual. Two densities of particle spacing were used in an effort to preclude cells 
voiding of particles in the leading edge flow. The densities were 36 and 100 per square 
cm. The interface shape for a tank with c = 0. 68, Bo = 10, and a fill of 25 percent was 
taken from Concus where the interface shape is graphically presented {Ref. 12). 

Marker particle plots for this case appear in Figure 13 and the vector plots in Figure 
14. The grid for this model tank is shown in Figure 14A. The cells range in size 
from 1 X 1 cm to 0. 125 x 0. 125 cm. Four of the smaller cells in a block make up the 
slit which represents the holes in the thrust barrel. This is discussed in greater 
detail in simulating the full-scale thrust barrel. The thrust barrel wall is made up 
of two rows of cells 0. 125 cm in thickness which are boundary cells and do not 
communicate information across the wall. In reviewing the make-up of the grid, 
it is now apparent that a finer subdivision of cells is needed outside the thrust barrel 
in the vicinity of the lower curved boundary. The single row of cells here which 
become filled with liquid did not provide adequate modeling since they were only one 
row in thickness and were all surface cells with empty cells above. However, when 
these were subdivided for Case 3, other problems in communication arose which will 
be discussed in Case 3. The inadequacy of arbitrary boundary modeling with thin 
films was uncovered in this case and was not solvable here or in Case 3. 

Flow progressed down the outer side wall in a satisfactory manner. The use of 
surface tension/contact angle logic on curved boundaries was demonstrated for the 
first time in a model simulation. Leading edge velocities were comparable to 
expected velocities, in a range of 0.8 to 1.0 of free-fall velocity. At 0.40 seconds, 
the maximum velocity vector was 20 cin/sec. No comparison with test data is 
available since the resolution of flow in the curved wall drop tower test containers 
is not of adequate clarity to present here; additionally, flow into the thrust barrel 
was not adequately defined in drop tower tests for quantitative comparison. After 
some difficulties were resolved with velocities on open cell faces near the curved 
boundary and the thrust barrel, flow did cover the lower hole and started to fill the 
thrust barrel. Flow moved up the side of the thrust barrel at a velocity of about 15 
cm/sec and passed above the thrust barrel at 0. 72 sec. The thrust barrel commenced 
filling at 0, 60 sec and continued to fill until the end of simulation at 0. 77 sec. 
Velocities through the slit were relatively constant during this period, ranging from 
12 to 16 cm/sec. At this rate of filling, the thrust barrel would be half full in 0. 5 
seconds from start of filling. This compares to a full-scale fill-time for half full of 
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Figure 14. Velocity Vector Plots for Model Centaur Oxidizer Tank Settling, Case 2 








approximately 25 seconds. The results of this simulation indicated the potential for 
Centaur oxidizer tank modeling, and simulation activities proceeded to Case 3, the 
full-scale tank. 

3.3 CENTAUR LO2 TANK SETTLING FULL-SCALE - CASE 3 

A new grid structure was developed to model the full-scale Centaur oxidizer tank. Fore- 
most, improved modeling of the flow along the lower bulkhead was required. Resolution 
cell-size was inadequate in Case 2. It was apparent that more cells had to be full, i. e. , 
the film had to occupy more than one layer. It was also desirable to modify the location 
of the thrust barrel side-hole in relation to the wall such that flow would not occur 
directly into the hole. It was correctly anticipated that the film along the wall in Case 3 
would be thinner due to lower initial liquid residual and that it would turn the comer on 
reaching the thrust barrel. It was also desirable to improve the modeling of the thrust 
barrel in line with moving the hole up the side. The dimensions of the thrust barrel are 
shown in Figure 12B. Three slits (annuli) were used to simulate the hole areas. 
moving them radially, the correct area could be achieved for top holes. The marker 
particle plots for this simulation are presented in Figure 15 and the velocity vectors in 
Figure 16. The full-scale grid configuration is shown in Figure 16A. The full-scale 
tank was shown earlier in Figure 11. 

The full-scale tank was represented with 390 cells and 430 particles. Cells ranged in 
size from lx 1cm to 0. 0625 x 0. 0625 cm, the latter used to represent the thrust barrel 
holes. A scale factor of 21. 7584 was used to scale the problem to full size. Top slits 
were made up of four small cells (0. 125 cm total) whereas the side slit was 0. 0625 cm 
(unsealed) and only a one by two horizontal array. As evidenced in Figure 16A, the grid 
in the vicinity of the lower bulkhead is essentially rectangular with good resolution of cells. 

The interface shape for the 9 percent liquid present with Bo = 15 was obtained from Concus 
(Ref. 12) by interpolation for shapes between 12 percent liquid and no liquid and between 
Bond numbers of 10 and 30. The flow was Initiated with surface tension logic active and 
a wall contact angle of 20° , which was measured from a surface marker particle at the 
center of the nearest boundary cell (Section 2. 2. 5). The simulation advanced down the 
sidewall in an uneventful manner. Velocities continued to increase proportional to the 
settling time. After ten seconds, some anomolies arose with the th in film thickness. 

Cells one away from the boundary were intermittently voiding and filling as particles 
moved in and out of them. From this point in time, certain cells had to be continually 
monitored to maintain them full without neglecting the restriction of maintaining a thin 
film. This procedure was required along the entire arbitrary boundary to the thrust 
barrel. Further difficulties developed in the transmission of advance cell velocity and 
pressure information when the film was only one cell thick in many instances. Extran- 
eous velocities occurred on empty cell faces which were unrealistic and raised doubt in 
the accuracy of the modeling of thin films such as experienced here. It appeared contin- 
uity in flow down the arbitrary boundary could not be obtained, an apparent shortfall in 
the code. After considerable attempts, solutions were obtained for flow into the comer 
area formed by the curved bulkhead and the thrust barrel. At 20 seconds, flow reached 
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Figure 15, Marker Particle Plots for Full-Scale Centaur Oxidizer Tank Settling, Case 3 
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16.0 SEC F.T= 18,0 SEC G, T = 19.4 SEC H. T = 21.0 SEC 

Figure 16. Velocity Vector Plots for Full-Scale Centaur Oxidizer Tank Settling, Case 3 









the thrust barrel and turned upward with a full-scale velocity of 20 cni/sec. It is reason- 
able to assume, from the results of Case 2, that thrust barrel filling would have occurred 
from this time onward. In light of the discontinuities experienced in the pressure and 
velocity fields along the boundary, the validity of the results for the pressure field 
causing flow through the hold are in doubt. The simulation was discontinued since 
computer runs would only proceed a few steps, less than one-half second, before some 
velocities along the lower boundary would exceed reasonable values. 

It is concluded that without further program modifications, the modeling of flow of 
extremely thin films of liquids over arbitrary boundaries should not be attempted. 

Further examination of program logic in these areas with respect to velocity and pressure 
field propogation is required. 

3.4 CENTAUR LH 2 TANK MIXING - FULL-SCALE - CASE 4 

A requirement exists during extended space missions to control cryogenic tankage pres- 
sure rise rates caused by fluid stratification. Most concepts currently under consideration 
utilize a mixer which consists of a jet pump to ensure the fluid in the tank is well-mixed. 

In reduced- gravity the mixer power requirement to pump liquid from one end of the tank 
to the interface is much smaller than in one-g. The strength of the jet does not have to 
overcome the pressure drop resulting from a gravity-produced liquid head. Mixing is 
also influenced in one-g by density differences, whereas in reduced-gravity these effects 
are minimized or absent. It is difficult, if not impossible, to maintain a similitude of 
the appropriate force and flow field variables for the two cases. Hence, one-g simulation 
of low-power mixing is not necessarily reliable. 

Only recently, with drop tower work at the NASA/Lewis Research Center, has low-g 
fluid jet mixing been Investigated experimentally. The marker-and-cell code provides 
a simulation tool to examine the full-scale mixing problem in reduced gravity fields. 

This simulation. Case 4, is representative of the potential for acquiring flow field data 
from such a simulation. The flow field is an adequate indicator of destratification. 

3.4. 1 SYSTEM SIZING - Previous test results in one-g with cryogenics and water to 
investigate de stratification for space missions afford one source of data for sizing the 
jet flow. The correlations for mixing from the chemical industry afford a second 
source of data. In a recent Convair effort (Ref. 13) the above sources were compared 
and the correlation which best represented available data was selected. To examine 
the variables and to show the basis of a design criteria, the work of these investigators 
will be briefly reviewed in chronological order. 

Fossett and Prosser (Ref. 14) studied the mixing of a Na2C02 solution in a tank with a 
jet mixer. The emf response of a pair of electrodes in the fluid to be mixed relative to 
a second pair in a homogeneous fluid were the measure of mixing, hi their data 
analysis, they determined a mixing factor which was independent of the jet Reynolds 
number. Their correlation for the mixing factor is a constant of 9 resulting in 
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(32) 


Dj /Dt^ = 9 

where 

Sjj^ -- mixing time in sec 

Vj = exit jet velocity, ft/sec 

Dj "jet exit diameter, ft 

Dt = tank diameter, ft 

The data from tests performed with four different configurations and fluids to analyze 
fluid mixing in aerospace applications is plotted in Figure 17 along with the above 
correlation. The correlation is adequate in view of data scatter but points up a 
potential inadequacy of their defined mixing factor. 

Recognizing this deficiency in the above correlation factor. Fox and Gex (Ref. 15) 
proposed a more complex factor to correlate their test results. Their experimental 
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Figure 17. Data Correlation for Selection of Mixing Factor, Equation by Fossett 
and Prosser 


38 




data was obtained from jet mixing of a known acidic solution into a tank of basic 
solution containing an acid-base indicator. Their correlation is presented in Figure 
18 where it is noted the laminar jet regime ends at Rej = 2000 and the turbulent region 
employs a different correlation. Only the turbulent regime is of interest here. Their 
work, in agreement with the above, suggests the parameter VjDj to be the key variable 
For the turbulent regime, they also include the liquid height and jet Reynolds number 


where 


e 

m 


(Vj g^/^ 

yi/2 


= 118 Re, 


- 1/6 


(33) 


g = the gravitational acceleration, ft/sec^ 

Y = liquid depth, ft 


The use of g as a low-g variable acceleration field rather than a conversion constant is 
questionable. Results in Figure 18 indicate a higher constant than 118 is required to fit 
the data, i. e. , mixing times are two-to-four-fold that predicted with this correlation. 



Figure 18. Data Correlation for Selection of Mixing Factor, Equation by 
Fox and Gex 
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The preferred correlation for mixing time is that presented by Okita and Oyama (Ref. 
16). Their mixing tests were also of a chemical nature and required establishing equal 
potential readings in NaCl solutions. In correlating results they defined a mixing 
factor which was independent of jet Reynolds number for Re > 5 x lO"^ of the form 


9 in Q 


D, 


1/2 


Vt d. 


= 5.2 


(34) 


where 


J 


= volumetric flowrate, Vj / 4 

2 


Q 

Vl = liquid volume ~ Y ff / 4 
Formulated in a manner similar to the above Investigators 
0 . 


^m (D, V ) 

77^. 


= 5.2 


(35) 


The data correlation is presented in Figure 19. Since It significantly reduces scatter and 
is the most conservative of the three correlations at high jet Reynolds numbers, it was 
used for defining mixer variables in the marker-and-cell Case 4 simulation. 
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Figure 19. Data Correlation for Selection of Mixing Factor, Equation by Okita and Oyama 
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REP RODUC: BDiTY OF TBK 




Another investigator evaluated the effect of the jet on the interface as it affects inter- 
facial mixing. Jet mixing was evaluated by Potli (Ref. 17) with water heated by wall 
heaters to develop a stratified layer which could be broken-up with a given strength 
jet. A criterion was theoretically established defining a critical Weber number such 
that mixer exit velocities in excess of this value would result in interface break-up. 

The correlation which was theoretically arrived at is 

Wecrit= V.^ D. //9= 1.055 Y / Dj (36) 

J 3 

for Bond numbers in the range -10 < Boj (. < 10. For a liquid height Y of 300 cm and 
Dj - 12 cm in the current study, the value for = 26 is independent of Bond 

number. 

There are, however, further considerations on the acceleration environment in the 
determination of ullage break-up from a submerged liquid jet. Poth comments that 
above a Bond number of 10, a Froude number criteria must also be considered (Ref. 

17). This Froude number criteria was derived for a tank with L/Dt of 2 for conditions 
for a submerged jet to reach the upper bulkhead. The criteria is presented (Ref. 17) 
in terms of Froude number 

^^crit = / g Dt = 0. 422 (D/Dj)^ (37) 

Considering the dimensions of the Centaur tank and the 12 cm jet outlet proposed above, 
the critical Froude number is 272. 

The critical Froude number criteria can also be expressed in terms of the jet critical 
Weber number as a function of the Bond number. With appropriate substitution and a 
change from tank diameter to tank radius for the Bond number (g 

WCcrit = 1- 688 (Dj./Dj) BOj,^ j. (38) 

The change in limits for the two correlations is Boj., t = 2. 5 rather than Bo = 10 above 
based on tank diameter. To consider this transition between the two correlations, they 
are plotted in Figure 20. It is apparent that a transition zone is required between the 
two correlations; however, the influence of the acceleration field is significant to a Bond 
number below 1. 

Some preliminary results from recent drop tower tests at NASA Lewis Research Center 
(Ref. 18) represent the only available low-g mixer data. Their work was performed in 
a 10 cm diameter cylinder with flow into the cylinder through an 0. 4 cm diameter jet 
colored by dye. Their results in very low-g verified the Weber number correlation for 
interface break-up. In tests where jet velocities did not exceed the critical Weber 
number, the flow patterns reached the interface and flowed laterally to the side wall and 
did not break the interface. The length of the drop did not provide adequate viewing 
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Figure 20. Weber and Froude Number Criteria for Interface 
Breakup 


time for dye trace evaluation of vortices or flow contact with the wall. Future tests 
at low thrust (finite drag Bond numbers) would provide the much needed verification of 
the critical Froude number criteria. 

In sizing the mixer for this full-scale Centaur fuel tank simulation, the criteria of 
Okita for mixing time and of Poth for a critical Weber number defining interface 
break-up were considered. The design procedure that follows is for a Bond number 
less than one using Equation 36 since the authors were unaware of the more restrictive 
criteria imposed at higher Bond numbers by the Froude number criteria. Since liquid 
depth above the mixer was specified as 300 cm with a tank diameter of 304. 8 cm, the 
two equations can be solved simultaneously if a factor a is introduced for the number 
of times above We^rit that is selected for operation, hi this series of trial and error 
calculations on two equations, the following variables are optimized: maximize 
turnover, minimize mixing time, maximize or to assure interface break-up within the 
constraints of a low velocity jet with a reasonable diameter, say 12 cm. It was 
preferred to simulate low velocity jets since high velocities reduce time-steps in the 
marker-and-cell technique. The chosen operational variables are summarized in 
Table II. 
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TABLE n. DESIGN VARIABLES FOR CENTAUR FUEL TANK MIXER - 
CASE 4 


Vl 

= 21. 6 X 10^ cm^ (60% liq.) 

Rej = 

1, 06 X 10^ 

Dt 

= 304,8 cm 

Wecrit = 

26.4 

Y 

= 300 cm 

Weact = 

132. 7* 

0 

= 27. 07 cmS/sec^ 

^^crit ” 

272 

Vj 

g 

= 12 cm 

F^act = 

56* 

= 17.3 cm/sec 

Q 

1956. 6 cm^/sec 

= + 0. 0174* cm/sec^ 

9m = 

2281 sec 

Bo 

= 15* 

®t 

11040 sec per volumetric 
turnover 

* It will be pointed out later that an 
finally used. 

acceleration field of 0. 116 with Bo of 100 was 


The data on volume turnovers is presented since this is a correlation variable for 
some authors. For a given diameter jet the mixing times of the first three authors 
above reduce to fractional volume turnover equivalents of 0.358, 0.073, and 0.207, 
respectively; i. e. this volume turnover has occurred in the predicted time for mixing, 
%n. The considerable extent of disagreement between workers is obvious. In 
Reference 13, Convair experienced tank homogeniety with LO 2 in a range from 0. 024 
to 0. 26 turnovers suggesting the non-reliability of this statistic. For the Case 4 
simulation this would suggest mixing times of 265 to 2870 sec. 

Once the mixer variables were determined, the details of establishing the grid for the 
problem were defined next. As indicated earlier, high velocities in a cell result in 
small time-steps because of time-step criteria in the computer code which limit 
particle travel to one-half a cell dimension. If the exit of the mixer, 6 cm radius, 
were used as a cell dimension with, for example, a 12 cm height, the time step 
would still be limited to 12 cm/(2 • 17.3 cm/sec) or 0.347 seconds. A solution to 
this was to use a large cell dimension for the mixer and, using the continuity equation 
and existing data on jet spread angle and entrainment, to determine the equivalent 
cell inlet and outlet velocities. From previous experience with the Centaur hydrogen 
tank, a cell-size of 152.4/24 cm in width (10, 886 cm) by 152.4/7 cm in height (21.772 
cm) was selected. 

The literature on axisymmetric submerged jets in an incompressible fluid provided 
data on the spread of the jet and the entrained fluid (Ref. 19). The jet spread angle 
is 15° . Thus within one cell length the 6 cm outlet radii can be represented by flow 
out a cell width. The fluid entrained by the jet is a function of the distance from the 
exit. These regions are defined as initial, which exists for 9 jet radii, a transition 
region, and a main region which starts at about 20 jet radii. In the initial region, 
the entrainment is determined from 
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( 39 ) 


_ _2 -4 

Q’ = (1 + 0. 073 X + 0. 002 x = . 0001 x + . 000002 x ) 

where 

Q' = volumetric flowrate at a distance x 

Qq = exit volumetric flowrate 

X = x/r^ with r^ nozzle radius 

x = centerline distance from the nozzle exit 

In the main region, the expression is simply 

Q = Qo (0. 155 X) (40) 

The first equation was used to determine entrainment within the mixer cell with Q'/Qq ~ 
1. 28. The equivalent inlet velocity at the top of the cell (Figure 21) is -5. 256 cm/sec, 
the cell side velocity for the cylindrical sidewall surface is -0. 377 cm/sec representing 
entrainment producing an outlet flow of -6. 763 em/sec. This latter value permitted a 
time step of 21. 77 cm/(2 • 6. 763 cm/sec) or 1. 61 sec. This was extended to 2 second 
steps since the flow in this region was predictable and cell voiding below the mixer was 
not a problem. 

3.4.2 MIXING SIMULATION - The physical characteristics of this full-scale LH 2 tank 
simulation were a 60% liquid-fill case represented by the grid configuration shown in 
Figure 23A and made up of 360 cells with initially 2407 particles. The basic grid 
dimension was 1 cm x 1 cm with a scale factor of 21.7714, with sub-divisions of this 

basic grid. The gravity level was initially -0. 0174 cm/sec2 
representing a Bond number due to drag of 15; however, this 
proved unsatisfactory and an acceleration of 0. 116 cm/sec 
(1. 2 X 10"'^ g’s) was selected "(Bo = 100). The problem of 
j non-convergence arose at the low Bond number and it was 
decided that the demonstration of mixing flow fields would 
also be of interest at a higher Bond number. This non- 
convergence problem exists due to the lack of well-defined 
' pressure field; comments on Case 1 apply here also. 

For the initial selected Bond number of 15, the critical 
Weber number is 640, and for the Bond number of 100 
which was used, Wej^crit ~ 4300. Since the Froude number 
is a significant factor in the ullage break-i^), the jet 
momentum was insufficient to break the interface. None- 
Figure 21. Simulated theless, mixing in low earth orbit remains to be of 

Mixer Cell considerable interest where drag Bond numbers in the 
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Figure 22. Marker Particle Plots for Full-Scale Centaur Fuel Tank 
Mixing Simulation in Low-g, Case 4 
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. T = 0.0 SEC B. T= 100 SEC C. T = 200 SEC D. T = 300 SEC E. T = 350 SEC F. T = 400 SEC 



G. T = 450 SEC H. T = s6o SEC I. T = 550 SEC J. T = 600 SEC K. T = 650 SEC L. T - 700 SEC 

Figure 23. Velocity Vector Plots for Full-Scale Centaur Fuel Tank 
Mixing Simulation in Low-g, Case 4 


46 




range 5 to 100 are applicable. With the low mixer exit velocities selected here based on 
the Weber number criteria at low Bond number, interface breakup could not be verified 
because the Froude number criterion was unsatisfied. The flow conditions modeled can 
still indicate the potential for effective destruction of stratification at the interface due 
to the interface velocities. 

A second problem in the initial set-up was the interface for a Bond number of 100. This 
interface appeared to be unstable, even with surface tension logic in use. Initial surface 
velocities (instabilities) developed near the wall independent of mixer flow which would have 
masked later mixer flow at the interface. It is assumed a balance of forces did not exist. 
This was resolved with the adoption of a flat interface as illustrated in Figure 22A. 

The simulation proceeded smoothly with cells acquiring velocities as the mixing front 
moved forward. The Marker particle plots are shown in Figure 22 for representative 
times from 0 to 700 seconds mixing time, initially, particles were added each step 
to three locations near the mixer to preclude cell voiding. It became apparent that 
other locations later in the problem also had to have particle addition to avoid empty 
cells which were really full. The initial addition provides a useful pattern similar to 
dye or smoke to follow the leading velocity front. This figure shows an interface 
instability which developed late in the problem and resembled a surface wave. The 
interface velocities by 700 seconds are sufficiently strong to ensure interface mixing. 

The velocity vectors for the same time period are presented in Figure 23. This group of 
plots are unique in that velocities are comparable since the maximum velocity on all of 
the plots (6. 0 cm/sec at the center of the mixer cell) was held constant. The velocities 
throughout the tank increase with time and approach a steady-state condition. The veloc- 
ity vectors Indicate the vortices formed in the tank. A well-defined path down the center- 
line and across the interface, as well as the ejq>ansion of the jet, is indicated. The 
nearly horizontal flow below the mixer raises an area of interest. Since the velocities in 
most of the tank are small relative to the mixer inlet velocity, a change in scale was 
introduced for duplicate plots after 300 seconds. This magnifies all vectors by 15 and 
provides improved resolution for examination of Interface velocities. These vector plots 
are shown in Figure 24, 

The magnitude of velocity vectors are plotted at specific tank locations as presented in 
Figures 25 and 26. These locations are considered to evaluate the approach to steady- 
state mixing flow conditions in the tank. The velocities along the centerline described in 
Figure 25 are well-behaved and remain relatively constant once flow is established at that 
location. The diminishing strength of the jet with distance from the nozzle is easily 
detected. In Figure 26, the velocities show larger perturbations due to interface instabil- 
ities, however a definite trend toward steady-state is observed. The absolute value of 
these vectors is not so important as the trend they present for the mixing capabilHy of 
this type of flow. An approximation of less than 450 seconds for the mixer flow to 
progress down the centerline and to the tank wall is a logical conclusion. A further 
conclusion is the adequacy of this flow to reduce stratification in a period less than 
predicted by the correlations of Okita. 
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Figure 24. Magnified Velocity Vector Plots for Full-Scale Fuel Tank 
Mixing Simulation in Low-g, Case 4 
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The marker-and-cell code was shown to be a useful design tool for the analysis of 
mixer location and mixer flow patterns. This simulation method affords an economical 
approach to the analysis of full-scale mixing flows. 
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4 

CONCLUSIONS AND RECOMMENDATIONS 


Extensive logic to consider surface tension forces which used a spline fit to a monotonic 
series of surface points was successfully added to the marker-and-cell code ERIE. This 
logic further includes end conditions which permit the simulation of contact angles on 
straight walls, as well as both concave and convex walls. The accuracy of this logic was 
verified for specific liquid interfaces, for test cases and in four simulation cases cover- 
ing a range of gravity levels. These four cases were the primary emphasis for this 
study. Case 1 involved the simulation of impulsive settling for a Centaur model LH2 
tank; Cases 2 and 3 respectively considered model and full-scale Centaur LO2 tank 
continuous settling; finally Case 4 simulated mixing in a full-scale Centaur LH2 tank. 

A detailed examination of convergence procedures in code ERIE was required to alleviate 
the extensive iterations which occur when the simulation does not contain a significant 
acceleration field and the pressure gradient in the tank is small. Two new convergence 
methods were concluded to be less effective than the existing method. Therefore, the 
current method was improved with the addition of an extrapolation technique for the 
pressure field. This was based on the magnitude of the change in three successive 
Iterates every 40 steps in the iteration. This procedure worked effectively in all cases 
at both high and low gravity levels although further improvements are still needed. 
Convergence at reduced gravity levels remains a major area for improvement in this 
code. It must be concluded that the solution to improved convergence which permitted 
the modeling of impulsive-g was not of imiversal application to cover acceleration 
environments required for mixing simulations. 

The impulsive settling Case 1 was a partial success in that simulation continued for 1.88 
seconds through more than half of the settling flow. An unexpected stall of the flow along 
the lower bulkhead for one-half second could not be explained. This delay and converg- 
ence difficulties prevented computation of this case to the desired 3 seconds. At that 
time, flow had stagnated and surface tension forces were dominant in returning the fluid 
to an equilibrium configuration. Our results did indicate an expected decay of velocity 
after the decrease in the acceleration field. The marker-and-cell code is ideally suited 
for this analysis in that the acceleration field can be easily changed. The new version 
of the code models the significant force of surface tension, and the flow simulation is not 
only graphically given, but the magnitude of local velocities may be assessed. 

The simulation of the Centaur LO2 tank indicated the code ERIE could be used for 
detailed simulation. The flow into the thrust barrel through ports represented a new 
dimension in modeling. The case made extensive use of the variable grid capability 
with a wide range of cell sizes, this case also made extensive use of arbitrary 
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boundaries. In Case 2 Centaur LO 2 tank drop tower model simulation was achieved 
through 0, 77 seconds which was adequate to demonstrate settling flow down the aft 
bulkhead and up the side of the thrust barrel. Flow into the thrust barrel through 
side ports was adequately simulated and an estimate of fill rates and times was 
determined. Flow difficulties along the aft curved boundary precluded the investigation 
of flow on top of the thrust barrel. 

In Case 3, the full-scale Centaur LO 2 tank simulation, flow down the aft bulkhead was 
achieved for a simulation of 9 percent liquid present for 21 sec when flow reached the 
thrust barrel. Flow simulation only reached the bottom holes of the thrust barrel. 

The thin film occurring with the little liquid present resulted in more complex 
simulation problems than with Case 2 for 30 percent liquid. Unrealistically high 
velocities were being defined on the aft bulkhead due to the thin film and continuation 
of the case was impractical. It was concluded that the logic in the program depends 
upon a communication of pressure and velocity through full cells which cannot be 
satisfied by a single row of cells containing fluid along an arbitrary boundary. 

Modeling was not improved with the use of finer cell divisions after Case 2. It is 
concluded that code modifications will be required before satisfactory modeling of 
thin films on arbitrary curved boundaries can be achieved. 

The results of Case 4 confirmed mixing e^ectations. The case demonstrated the 
capability of the marker-and-cell code ERIE to simulate reduced-gravity mixing. The 
simulation was stopped at 700 seconds when steady-state conditions were realized. 

The flow patterns developed in a logical manner and reached steady-state in a time 
less than was predicted from available correlations for total mixing. The velocity 
vectors along the interface were strong and indicated the mixer flow had reached the 
wall prior to 700 seconds. This flow from the mixer to the interface would result in 
destratification. Theoretical predictions for interface break-up indicated jet 
momentum would be insufficient to break-up the interface and this was proven to be 
true. Some minor wave motion and Instabilities approaching the wall were character- 
istics of the final interface behavior. The Inability to gain convergence on mixing at 
low Bond numbers is not understood with respect to the results of Case 1. The mixing 
correlations appear to be conservative if the criterion is pressure decay caused by 
cold fluid from the mixer replacing the interface fluid. Earlier investigations were 
concerned with total tank homogenity. Case 4 did demonstrate the capability of this 
code to model mixing at a low cost. 

It is recommended that additional work be done to improve the convergence of this code 
for low-g simulation where the pressure gradient is small. This would be productive 
in examination of mixing at lower Bond numbers and with stronger jets; i. e. , higher 
jet exit velocities. The final Case 4 would not converge at low gravity levels. There- 
fore, a higher gravity level was attempted for Case 4 than in the impulsive-g Case 1. 
The temporary solution was to increase the gravity level from Bo = 15 to 100, however, 
this would not be satisfactory for cases which should be modeled at Bond numbers 0. 1 
to 10. The results of impulsive-g settling indicate the promise of the efficiency 
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of such a settling approach. This code can examine sequencing and optimum thrust 
application for full-scale simulations. The new capability of the code with surface 
tension forces opens several areas of low-g fluid motion for examination: impulsive 
settling, fluid mixing during coast and surface tension driven flows. 
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APPENDIX A 


NOMENCLATURE 


Bo 

c 

D 

D 

EPS 

G 

Sr 

Sz 

h 

n 

P 

Q 

r 

R 

^ 1 ’ ^2 

Re. 

3 

t 

u 

u 

V 

V 

Vm 

V 


Wej 

X 


Bond number, g^, rj. / ^ 
wave speed 

velocity divergence = (1/r®) (3 r®u/dr) + (dv/3z), distance 
diameter 

convergence criteria 
cell flag 

radial acceleration 
axial acceleration 
height of liquid 

unit normal defining a boundary segment 
pressure, L^T“^ 
jet volumetric flow rate 
radial coordinate 

radial distance to center or side of cell determined by subscript 
radius of curvature 
jet Reynolds number 
time 

radial component of velocity 

radial storage variable, radial component of the tilde velocity 
axial component of velocity 

axial storage variable, axial component of the tilde velocity 

liquid velocity at midpoint of boundary segment computed with 
area weighting interpolation scheme 

velocity, volume 

jet Weber number, D.//S 

horizontal coordinate in rectangular coordinate system 




■V* NOT FILMED 



y vertical coordinate in rectangular coordinate system 

Y liquid height 

z axial coordinate in cylindrical coordinate system 

a geometric parameter, a- 1. 0 in cylindrical coordinates and equals 

0. 0 in plane (cartesian) coordinates 

6r incremental step in the r direction 

6 z incremental step in the z direction 

8 kinematic surface tension, T"^ 

€ boundary sensing parameter 

9^ mixing time, T 

X minimum mesh dimension, minimum of 6 r or 5 z 

u dynamic viscosity 

T) ratio of lengths in defining surface pressures 

u kinematic viscosity, 

P density 

o surface tension 

5 relaxation parameter 

0 true pressure normalized to unit density 

if arbitrary pressure normalized to unit density (pseudopressure), 

u) vortieity 

Superscripts 

k iteration index 

n time cycle index 

Subscripts 

B bottom 

L liquid 

j iet 

N identifier for primary cell under consideration 

NL left boundary of cell N 

NB bottom boundary of cell N 

T topi 

t tank 
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